Diffusion-limited reactions on a two-dimensional lattice with binary disorder 
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Reaction-diffusion systems where transition rates exhibit quenched disorder are common in physical and 
chemical systems. We study pair reactions on a periodic two-dimensional lattice, including continuous deposi- 
tion and spontaneous desorption of particles. Hopping and desorption are taken to be thermally activated pro- 
cesses. The activation energies are drawn from a binary distribution of well depths, corresponding to 'shallow' 
and 'deep' sites. This is the simplest non-trivial distribution, which we use to examine and explain fundamental 
features of the system. We simulate the system using kinetic Monte Carlo methods and provide a thorough un- 
derstanding of our findings. We show that the combination of shallow and deep sites broadens the temperature 
window in which the reaction is efficient, compared to either homogeneous system. We also examine the role 
of spatial correlations, including systems where one type of site is arranged in a cluster or a sublattice. Finally, 
we show that a simple rate equation model reproduces simulation results with very good accuracy. 

PACS numbers: 98.38.Bn, 68.43.-h, 98.38.Cp 



I. INTRODUCTION 

Reaction-diffusion systems are successful models to de- 
scribe a large variety of phenomena in physics, chemistry, and 
biology [1,2]. They may involve one or more reactant species 
that diffuse and react with each other on a surface or in the 
bulk. In particular, surfaces often catalyze chemical reactions 
between adsorbed atoms and molecules. The densities of the 
adsorbed chemical species and their reaction rates depend on 
parameters of the surface and on the temperature. Microscop- 
ically, one can describe the diffusion of particles on the sur- 
face as a random walk between adsorption sites. In homoge- 
neous systems all the adsorption sites are identical. However, 
most systems are heterogeneous, involving different types of 
adsorption sites with a broad distribution of binding energies. 

The present study is motivated by a specific example of an 
important surface process, namely the formation of molecular 
hydrogen on dust grains in the interstellar medium [3-8]. Hy- 
drogen atoms impinge from the gas phase onto a grain, and 
diffuse on its surface. They may either desorb thermally from 
the surface, or encounter each other and form a molecule. This 
defines a reaction-diffusion system in a spatially confined re- 
gion. In this article we are concerned with steady-state sys- 
tems, when the hydrogen recombination efficiency is defined 
as the fraction of impinging particles that end up (and eventu- 
ally desorb) in molecular form. This efficiency plays an im- 
portant role in the evolution of interstellar clouds. Typically, 
there is a narrow window of temperatures in which recombi- 
nation is efficient. At lower temperatures, the atoms are not 
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sufficiently mobile to react, whereas at higher temperatures 
they desorb too quickly. 

Assuming that all rates are spatially homogeneous, the sys- 
tem is well understood analytically. A zero-dimensional mas- 
ter equation for the particle number distribution [9, 10], to- 
gether with a proper definition and calculation of the reaction 
rate coefficient in terms of a first-passage problem [11, 12], 
suffices to accurately describe the many -particle system [13]. 
It is very important, however, to consider disorder in the local 
rates of hopping and desorption of the particles. As we alluded 
to earlier, this is not only of theoretical interest. In fact, in the 
astrophysical context, the disordered case is much more real- 
istic, and it is long known that disorder potentially enhances 
the efficiency dramatically [5]. However, the combination of a 
confined two-dimensional region, rate disorder and the many- 
particle reaction-diffusion dynamics makes this problem no- 
toriously hard to tackle analytically. Kinetic Monte Carlo 
(KMC) methods can be used to simulate such systems [e.g. 
14, 15], and algorithms are still subject to improvement [16, 
which also compares related approaches]. They remain com- 
putationally expensive, however, and a systematic understand- 
ing of the effects of disorder is still missing. 

Here we start such an analysis for the simplest form of rate 
disorder, where each lattice position corresponds to either a 
standard ('shallow') site, or to a strong-binding ('deep') site, 
with enhanced binding energy. While we strive to keep this 
a theoretical self-contained work, our models and questions 
are motivated by applications and should easily translate to 
practice. This is one reason why we have chosen thermally 
activated rates throughout, and present most results in terms 
of temperature, and on scales relevant to the astrophysical 
problem just described. In the latter context, our work is rel- 
evant to systems combining physisorption (shallow sites) and 
chemisorption (deep sites) [ 17] , aside from features particular 
to specific material systems. Using such a discrete distribution 
turns out to be conceptually different from the case of con- 
tinuous distributions of binding strength [14], in which well 
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depths drawn from tails of the distribution may significantly 
affect the temperature window of high efficiency. 

Our goal in this paper is to provide a thorough under- 
standing of all relevant mechanisms of the described reaction- 
diffusion system. Most importantly, if we start from homo- 
geneous systems of either standard or deep sites, their tem- 
perature windows of high efficiency will typically be sepa- 
rated by a gap. It is a natural question whether a mixture of 
the two types of sites still exhibits two separated peaks, or 
whether (and under what conditions) the efficiency is high for 
in-between temperatures. 

Our findings are relevant for other systems as diverse as 
catalysts [18], exciton trapping in photosynthesis [19], exciton 
transport in semiconducting nanosystems [20], and diffusion- 
limited reactions on biomembranes [21]. The generalization 
of our results to these and other related contexts should be 
straightforward. 

The paper is organized as follows. In Sec. II we define the 
system and our notation. The following Sec. Ill provides a 
qualitative picture which identifies three temperature regimes 
and describes the relevant processes in each. In Sec. IV we 
give a systematic account of extensive KMC simulations and 
discuss the observed behavior in detail. This includes the 
study of spatial correlations in the quenched disorder. Sec- 
tion V presents a simple yet accurate rate equation model, and 
we explain the difference to the homogeneous case. We derive 
an expression for the efficiency in the most interesting regime. 
Finally, we present our conclusions in Sec. VI. 
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FIG. 1: One-dimensional cut through the energy landscape of our 
model. 



way to achieve this is by choosing W\/ai = W2/C12, or equiva- 
lent^, £Vj —E a] = Ew 2 —E U2 , and we will employ this choice 
throughout. The number of sites visited by a single particle 
before desorption becomes then independent of disorder. 

To establish a connection to surface chemistry problems, 
we think of type-1 sites as standard or 'shallow' adsorption 
sites, and of type-2 sites as strong-binding or 'deep' sites, with 
Ew 1 > Ew { ■ A one-dimensional cut through such an energy 
landscape is sketched in Fig. 1 . 



III. QUALITATIVE DISCUSSION 
A. Homogeneous systems 



II. MODEL AND DEFINITIONS 

We consider a system of a single particle species on a two- 
dimensional square lattice of S sites with periodic boundary 
conditions. Each lattice site is characterized by a binding en- 
ergy, which can take one of two values — we call this a bi- 
nary lattice. The number of sites of either type is denoted 
by Si (i = 1, 2), and S = Si + 5V Particles impinge onto the 
lattice at a homogeneous rate / per site. If a site is already 
occupied, the impinging particle is rejected. In the context 
of surface chemistry this is known as Langmuir-Hinshelwood 
(LH) rejection [22]. 

Particles explore the lattice by hopping to neighboring sites 
with an (undirected) rate a, and they can desorb from a site 
with rate W. Both rates depend on the binding energy at the 
particle position. If two particles meet on one site, they form 
a dimer and leave the system immediately. The key quantity 
of such a system is the efficiency 77, defined as the ratio be- 
tween the number of particles that react and the total number 
of impinging particles, when the system is in a steady state. 

In view of possible applications, we choose rates to be ther- 
mally activated by a system temperature T. The activation 
energy for desorption is denoted Ew r which we identify with 
the binding energy at the particle position. Similarly, hopping 
from a type-/ site has an activation energy E clj . All rates share 
the attempt frequency v, so that, e.g., Wj = vexp (—EwjT) — 
here and in the following energies are measured in tempera- 
ture units. We want to ensure detailed balance. The simplest 



For a homogeneous lattice, the dependence of the efficiency 
on the temperature, t]{T), is known [9, 10]. At low tempera- 
tures the particles are nearly immobile, thus they do not meet 
other particles. Therefore, the lattice is highly occupied, in- 
coming particles are mostly rejected, and the efficiency is low. 
For higher temperatures, hopping processes are activated and 
the particles begin to explore the lattice. This leads to more 
frequent encounters, so the efficiency rises. When the tem- 
perature is increased even further, the particles tend to desorb 
before encountering each other, the lattice coverage becomes 
small, and the efficiency decreases. 

For the homogeneous system, the corresponding tempera- 
ture bounds have been obtained using rate equations [10]: 
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is the temperature above which the kinetics becomes second- 
instead of first-) order, whence desorption ends the typical 
particle residence and the efficiency is low, and 
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is the temperature below which particles arrive faster than they 
hop, leading to dominant LH rejection and low efficiency. The 
average of these two bounds reads 
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and corresponds to the temperature of maximum efficiency. 

If the binding energy Ew is increased, the efficiency max- 
imum is shifted towards higher temperatures, and the shift is 
directly proportional to the change in binding energy. For a 
binding energy difference AE = Ew 2 — E^ l of two (otherwise 
equal) lattices of either type-1 or type-2 sites, the relation be- 
tween the temperatures of maximal efficiency is given by 
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whereas the peak width 7j up — 7J low is the same. 



B. Binary systems 

Now we consider the binary lattice introduced in Sec. II, 
with binding energies Ew l and Ew 2 ■ To each site, we randomly 
assign a binding energy. There is a typical length for a parti- 
cle to find a strong-binding site. This length obviously short- 
ens when there are more and more of these sites on the lat- 
tice. At low temperatures around the efficiency maximum of 
the type-1 sites, particles can only diffuse on and desorb from 
these shallow sites, while particles landing on or hopping onto 
strong-binding sites cannot leave by hopping or desorption, 
since the binding energy is too high. Recombinations either 
take place on the shallow sites, or by hopping to an occupied 
neighboring strong-binding site. For very high temperatures 
around the efficiency maximum of the deep wells, the particles 
diffuse on, desorb from and recombine on those, while on the 
shallow sites, they desorb too quickly to allow any other pro- 
cesses. But in the intermediate temperature regime — right 
of the shallow peak, left of the strong-binding peak — some- 
thing different happens. Here, the temperature is too low for 
dynamics on deep wells, so particles encountering a deep well 
are stuck. On the other hand, particles on shallow sites tend to 
desorb rather quickly, and thus do not recombine on such sites. 
But if they find a deep well before desorbing, they are trapped 
until another adatom shares their fate and they recombine. 

The simple random walk with traps has been studied exten- 
sively [e.g. 19, 23]. To leading order, the average number of 
steps a random walker performs before trapping is given by 
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where 52 is the number of deep wells and 5 is the total number 
of sites on the lattice. This leads to a trapping length 
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On the other hand, the typical radius of the area a walker ex- 
plores on standard sites before desorption is the random walk 
length [12] 
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Trapping now competes with desorption from shallow sites; 
the former only depends on the number of traps 52, while the 



latter is a function of temperature. As long as the random 
walk length 4w is larger than the trapping length £ tlap , the par- 
ticles are — on average — trapped before they can leave the 
lattice. For a given number of traps 52 this implies a high ef- 
ficiency approximately up to the temperature T eq where both 
lengths become equal. If this temperature lies above the inter- 
mediate temperature range where both pure systems have poor 
efficiency, we can expect a high efficiency throughout, hence 
a full 'bridging' of the gap. Since the efficiency is high over 
this whole temperature range then, we call this an efficiency 
plateau. We will calculate the value of the efficiency on such 
a plateau in a rate equation model in Sec. V C. 

Summing up, we can divide the temperature axis into three 
regions. The lowest temperatures where only particles on 
shallow sites are mobile, the intermediate regime where the 
particles behave like random walkers on a lattice with traps, 
and the high temperatures where particles become mobile on 
strong-binding sites. We now check this qualitative picture 
with KMC simulations. 



IV. KINETIC MONTE CARLO SIMULATIONS 
A. Setup 

In order to test our predictions, we carried out extensive 
kinetic Monte Carlo simulations. The standard algorithm pro- 
ceeds as follows [cf . 24, for a review] . We keep track of the 
full microscopic dynamics of continuous-time random walk- 
ers [25] with standard exponential waiting time distributions. 
In each simulation step, the current system configuration de- 
termines the list of possible elementary processes and their 
rates. By comparing a random number with the normalized 
partial sums of these rates we find the process to execute next. 
The simulation time is then advanced according to the total 
sum of rates and the configuration is updated. 

For a given realization, we wait for the system to reach the 
steady state before we measure the efficiency over 10 6 im- 
pingements. We use a square lattice of 5 = 100 x 100 sites. 
We choose the other model parameters inspired by an exem- 
plary system in the astrophysical application, to show the rel- 
evance of our work in this field, and since the corresponding 
system is known to exhibit interesting kinetic regimes. The 
flux of hydrogen atoms per unit surface area depends on gas 
density and temperature. The flux per surface site is given 
by the ratio between the flux per unit area and the density 
of adsorption sites, hence it depends on the surface morphol- 
ogy. More precisely, the flux per site is given by f — pv/(4s), 
where p is the density of hydrogen atoms in the gas phase, 
v is their average thermal velocity and s is the density of ad- 
sorption sites on the surface. To obtain typical values we use 
p = 10 cm~ 3 , v = 1.45 x 10 5 cm/s (which corresponds to a 
gas temperature of 100 K) and s = 5x 10 13 cm~ 2 which is the 
measured density of adsorption sites on the amorphous car- 
bon sample studied in Ref. 26. This results in a flux per site of 
/ = 7.3 x 10~ 9 s _1 . For the attempt frequency we choose the 
standard value of 10 12 s _1 which is commonly used through- 
out surface science. With each site we associate either the 
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standard binding energy Ew t = 658 K, as found for hydro- 
gen atoms on amorphous carbon [27], or an enhanced energy 
E Wl = E Wl + AE with AE = 250, 750 or 1500 K. The activa- 
tion energy for hopping reads E a] =511 K or E ai = E ai + AE, 
respectively. 

In each case, we determine the efficiency as a function of 
the temperature T, as well as of the relative frequency of 
strong -binding sites S2/S. We do this for up to four differ- 
ent ways of distributing the binding strengths. For dynam- 
ics with nearest-neighbor hopping of the particles, we either 
randomly assign to each site a binding energy with proba- 
bilities p\ and p2 = 1— pi, respectively, or we arrange the 
strong-binding sites in a regular sublattice, or we concentrate 
all strong-binding sites in a single square cluster. In the case 
of random assignment, S2 is then binomially distributed with 
parameter p%. To eliminate the fluctuations in S%, we average 
the efficiency over 20 realizations. In the following discussion, 
we can therefore identify Sj/S with its average /?, . For compar- 
ison with the rate equation model to be introduced in Sec. V, 
we also implement another kind of dynamics ( 'longhop ' case), 
namely hopping from any site to any other site of the lattice. 
This switches off any spatial correlations between the lattice 
sites and thus is best suited for comparison with an effective 
zero-dimensional model. 

Binary disorder models very similar to random assignment 
and the clustered case have been simulated before [14]. The 
authors were predominantly concerned with showing that 
such models can exhibit efficient reaction over a broader range 
of temperatures than homogeneous systems. Here we extend 
these findings to a systematic picture for the effect of the deep- 
site fraction pi and the energy gap AE. More importantly, we 
provide detailed explanations and analytic results which ex- 
plain all notable features of the simulation outcome in terms 
of microscopic physical processes. 



B. Results 

Figure 2 shows the results of our simulations. For each AE, 
we simulated systems with 1, 4, 25 and 50% of strong -binding 
sites. The random distribution is probably the most interesting 
regarding applications. Following the series of Figs, for each 
AE, we observe that the intermediate temperature regime is 
bridged in each case. This is in accordance with the analytic 
prediction of Sec. IIIB, since already for moderate deep-site 
fraction, the trapping length £ tmv is smaller than the random 
walk length l^, for all intermediate temperatures. The obser- 
vation holds at least up to AE = 2500 K (not shown), which is 
the largest value of AE that we have considered; beyond this 
energy scale one enters the regime of chemisorption, which is 
not our focus in this work. The bigger the difference of the 
binding energies, the more strong-binding sites are needed to 
form a genuine plateau, where the efficiency does not depend 
on the temperature. This complies with the ideas of Sec. Ill; 
when the deep-site peak is shifted to higher temperatures, T eq 
has to increase to warrant formation of a plateau. This is 
achieved by increasing the deep-site fraction. The variance of 
the efficiency between different realizations of random land- 



scapes was found to be negligible throughout. We also exam- 
ined the longhop case on such landscapes, and found that the 
efficiency varies just as much. Since this cannot be affected 
by any spatial correlations, we conclude that this variation is 
always due to fluctuations in the number of deep sites 52 only. 

The arrangement of strong-binding sites in a sublattice per- 
forms slightly better, compared to random assignment. This 
is not astonishing since the sublattice optimizes the distance 
between the traps. In the random case, small clusters of strong- 
binding sites can occur, in which a single trap is less efficient. 
An alternative picture is that the capture zones of individual 
traps typically have an overlap, which is minimized in the sub- 
lattice case. 

On a lattice with a single square cluster of strong-binding 
sites, there is no bridging effect for any energy difference 
or frequency of strong-binding sites. For high frequencies 
of either shallow or strong-binding sites, only one restricted 
peak emerges, while for intermediate frequencies of strong 
and shallow sites two nearly separated peaks appear. The ef- 
ficiency does not drop to zero in the intermediate temperature 
regime, because an exchange between shallow and deep sites 
takes place along the boundary of the cluster. However, since 
the boundary length scales as y/S, the fraction of boundary 
sites decreases with increasing S, and correspondingly the sup- 
pression of the efficiency in this regime becomes even more 
pronounced for larger systems. We checked this for a system 
of 500 x 500 sites (not shown). This is in contrast to the well- 
mixed case, where a finite fraction of sites are boundary sites 
(see below). 

For the longhop case, we first verified that results on a sub- 
lattice and a cluster landscape coincide, ensuring the correct- 
ness of the algorithm. The efficiency for this kind of dynamics 
outperforms even the sublattice results for nearest-neighbor 
hopping. This is because in the sublattice case, there is still 
the necessity for a particle to actually travel to a trap instead 
of having a non-zero probability to reach a trap on every step. 
A further analysis of this model is provided in Sec. V. 

In addition to our qualitative explanations, we numerically 
examine the dependence of the plateau efficiency value on the 
number of deep wells. First we note that for our choice of 
parameters, the efficiency value at J™ 3 * « 14 K always corre- 
sponds to the plateau value. From the results for AE = 750 K 
shown in Fig. 3 we infer that the way of distributing the strong- 
binding sites is of crucial importance. In the case of a sin- 
gle square cluster of deep wells, the efficiency decreases lin- 
early as 1 —S2/S, while for the random distribution the ef- 
ficiency first decreases more slowly (for less than 50% of 
strong-binding sites) and faster to the end (more than 50%). 
We propose that this effect is related to the border length be- 
tween shallow and deep sites, and use this connection to de- 
rive an empirical formula for the plateau efficiency. For ran- 
domly distributed deep wells, we calculate the border length L 
as function of Sj/S (cf. Fig. 4). We find a shallow site next to 
a deep site with probability (S2/SXI — S2/S). Since the orien- 
tation of the pair does not matter, we gain an additional factor 
of 2. Furthermore we have 2S possibilities to place such a pair 
of sites on a square lattice with S sites and periodic boundary 
conditions. So we find the following expression for the border 



5 



Horn, systems 
"T" 




Hom. systems 



Horn, systems 



10 12 14 16 18 20 22 24 




10 12 14 16 18 20 22 24 

4% 




10 12 14 16 18 20 22 24 

25% 




10 12 14 16 18 20 22 24 

50% 




10 12 14 16 18 20 22 24 




10 15 20 25 30 



FIG. 2: (Color online) Efficiency versus temperature for various fractions of deep sites. Left column AE = 250 K, middle AE = 750 K, right 
AE = 1500 K. Randomly assigned energies (blue line, diamonds), longhop dynamics (red line, circles). Only for AE = 750 K: sublattice 
(orange dashed line, crosses), and cluster (black line, squares). Vertical green line at T eq . The first row shows the results for homogeneous 
systems of only standard or only deep sites, respectively. 



length between shallow and deep sites multiple of this border length yields 

A7]=CL, (9) 
L = 4S ' ~s ( ! ~ ~s) ' (8) with C= (1.487 ±0.019) x 10~ 5 or 

Fitting the efficiency difference A77 = 77 random - 77 cluster to a ?7random ~ ' V + ( °' 595 ± °- 008) f ) ( 10) 
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FIG. 3: (Color online) Efficiency as function of the deep-site frac- 
tion for T = 14 K and AE = 750 K, for clustered deep sites (black, 
squares, r) c i uster ) and randomly assigned energies (blue, diamonds, 
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FIG. 4: (Color online) Border length L (green dashed line, left axis) 
and efficiency difference Ar/ (dark red, right axis), as function of 
deep-site fraction, for T = 14 K and AE = 750 K. Vertical axis scal- 
ing taken from data fit. 



for the empirical plateau efficiency value. The quality of the 
fit Atj °< L for KMC results underlines the role of the border 
length, and this corroborates our picture that the dominant re- 
action process on the plateau is by hopping from standard to 
deep sites. Further insight into the origin of Eq. (10) will be 
provided below in Sec. V C. 



V. RATE EQUATION MODEL 

Rate equations have been used previously to study reactions 
on finite surfaces with different types of sites. Using surfaces 
with varying roughness, where binding energies at a site are 
given by a vertical bond strength plus an additional lateral 
bond strength per in-layer neighbor of the landscape, KMC 
simulations were performed [15]. A rate equation model was 
then used to check that such a rough landscape model is con- 
sistent with surfaces deemed astrophysically relevant and ex- 
amined in the laboratory. To this end, the rate equations with 



standard energy parameters were time-integrated to predict 
the results of TPD experiments. 

Here we apply a rate equation model to quantitatively repro- 
duce our KMC findings as well as to further our qualitative un- 
derstanding of the system's behavior. From the definition of 
Sec. II we derive a set of rate equations for the total number Ni 
of particles on sites of type ;. Terms for desorption and influx 
are easily written down, whereas the reaction terms are more 
subtle, partly since in general, the reaction is not an elemen- 
tary process with given rate. For the time being, we denote 
the appropriate rate coefficients as A,, and refer to Sec. V A 
for details. The rate equations then take the form [15, 28] 



Si =f(S 1 -Ni)-W { Ni -A l Ni(S 2 -N 2 )-A 1 NiN 2 



N, 



(11) 



-2A l Nf +A 2 N 2 {S l -Ni) -A 2 N\N 2 , 
- f{S 2 -N 2 )- W 2 N 2 - A 2 N 2 {Si -Ni)- A 2 N\N 2 
- 2A 2 Nj +A 1 N 1 (S 2 -N 2 )-A 1 N 1 N 2 . 



Here the first two contributions cater for the impingement flux 
with rejection and the desorption of particles. For clarity we 
separated the remaining terms. The next two terms describe 
leaving to a site of the opposite type (either to an empty or 
to an occupied site). Then we account for reactions inside 
one population due to hops between sites of the same type, re- 
moving two atoms. The remaining two contributions describe 
gaining a particle by a hop from the other site type, and finally, 
losing one particle due to the reaction with a particle coming 
from the other population. 

It is tempting to substitute the 'internal' reaction term 
2AjNf by 2AjNi(Nj — 1), since it should really depend on the 
number of pairs. This is not adequate: In the rate equation 
treatment the Ni are continuous and can drop below unity, such 
that the reaction term (which we are ultimately interested in) 
could then become negative. The assumption that the reaction 
rate can be written as above is at the heart of the rate equation 
approach ("mass action law"). Equations (11) are easily de- 
rived from the full master equation using this assumption in 
the forms (Ni(Ni - 1)} « Nf and (NiN 2 ) w NiN 2 (where the 
expectation is over the joint probability distribution P(Ni,N 2 ) 
and the r.h.s. Ni's are already the mean values as above). 

The reaction terms also provide the recombination rate of 
the process. Adding up all terms proportional to the A, in 
N/t = Nj/t + N 2 /t, mere hopping terms (not leading to a re- 
action) cancel. Using that the reaction consumes two particles, 
we obtain the rate at which particles are removed by the reac- 
tion as 



2R = 2A x Nl + 2A 2 Nl + 2{Ai +A 2 )N l N 2 , (12) 

which can be simplified to 2R = 2(A i N l +A 2 N 2 )(N l +N 2 ). 
Relating this to the particle influx f(S\ +S 2 ) = fS gives the 
efficiency rj=2R/(fS). 
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A. The reaction rate coefficient 

The homogeneous system was treated analytically by rate 
equations [10, 27, 29], the master equation [9, 10], and mo- 
ment equations [30, 31]. For these methods just as for stochas- 
tic or numerical methods based on these approaches, the re- 
action rate coefficient is a crucial quantity, typically approxi- 
mated as A m a/S [32]. We have argued elsewhere that this 
neglects the nature of two-dimensional diffusion ("back diffu- 
sion") as well as the fundamental first-passage problem, the 
competition between a meeting (hence reaction) of particles 
and the prior desorption of a reactant. Hence we put some 
effort into a proper definition and evaluation of A [11, 12], 
and we claimed that these results should be applied in all men- 
tioned frameworks, including the rate equation treatment [13]. 

Here, we return to the choice A; = at/S, since the situa- 
tion is different. In rate equations such as Eqs. (11), there is 
no way to genuinely incorporate any spatial structure. This 
holds true for all zero-dimensional approaches, e.g., the mas- 
ter equation as well. However, in the homogeneous systems 
studied before, this neglect only concerns the spatial correla- 
tions in the particle residence probability, with well-studied 
effects [10, 11, 33]. In the heterogeneous system with its sep- 
arated populations, this approach additionally neglects site 
type correlations. 

For consistency, we are then forced to assume that a particle 
can reach any other site by a single hop. In particular, it hops 
to a site of type ; with probability Sj/S, and it meets a particle 
on an /-site with probability Nj/S. The conventional choice 
A, w at/S thus arises naturally if we use rate equations to de- 
scribe a system with site disorder, and we adopt this choice in 
the following. For a system with quenched spatial structure 
and nearest-neighbor hops only, this description corresponds 
most closely to the well-mixed case. 

B. Comparison with KMC simulations 

The rate equations (11) are exactly solvable at steady state 
by finding the real positive root of a third-order polynomial. 
However, the results are cumbersome and less than illumi- 
nating. We therefore directly opted for a numerical solver 
throughout. 

We find that the rate equations for the binary system repro- 
duce the outcome of extensive KMC (longhop) simulations 
for a wide parameter range of practical relevance to excellent 
accuracy (see Fig. 5). As noted in prior work [1 1], however, 
since we present our results as functions of temperature and 
parameters are thermally activated, we typically have rather 
steep rises or declines, when even factors of two or three in 
the efficiency need not appear substantial. This hardly ex- 
plains the overall accuracy, especially on plateaus and mod- 
erate peaks for r\ considerably smaller than unity. 

Results on the validity of rate equations to describe the 
model in the homogeneous case have shown that confinement 
to a finite surface renders the discreteness of particles and 
fluctuations in the particle number important [10, 11, 33-36]. 
Consequently, the mean-field approach of rate equations con- 



siderably overestimates the recombination efficiency in small 
systems. We do not see such effects for several reasons. We 
are interested in the behavior of the system with a substantial 
number of particles, when the effects of discreteness and of 
fluctuations in this particle number are strongly reduced. Fur- 
ther, the confinement of particles to a finite surface is also 
far less important than for the homogeneous system, because 
the majority of these particles is trapped in deep wells in the 
regimes of most interest, anyway. Finally, our system can- 
not be considered small, and we cannot preclude completely 
that differences might be more pronounced for smaller system 
sizes or different activation energies. 



C. Plateau efficiency 

A key question of this work concerns the bridging between 
the two efficiency peaks corresponding to homogeneous sys- 
tems of one type of sites. We have found a convincing quali- 
tative picture before in Sec. Ill, and we observe an efficiency 
plateau between the two (virtual) peaks for a wide range of 
conditions. What is the value of the efficiency along this 
plateau? 

We need to further simplify our rate equation model to ar- 
rive at a simple analytic answer. Several efforts to derive 
this from first principles were not met with success, hence we 
start from some observations: Figure 5 shows that whenever 
a plateau emerges in the efficiency, practically all recombina- 
tions are due to hops between the two types of sites. Only 
for very low concentrations of deep wells and when the ef- 
ficiency is close to unity, the lower temperature end of the 
plateau also includes a substantial contribution from recom- 
bination on standard sites. On the high-temperature end, any 
sizable contribution from reactions on deep sites already re- 
sults in an efficiency peak atop the plateau value anyway. 

For the plateau efficiency, we can therefore capture the 
essence of the model accounting only for reactions between 
the two populations. Type 1 denotes the standard sites, so 
we clearly have A\ ^ A^. We use this to neglect all terms 
proportional to Aj_, since they are small compared to their A i 
counterparts, but keep all flux and desorption terms. Our rea- 
soning is to retain as many terms as possible, to remove those 
for recombination inside the Nj populations, and since we can 
neglect the reaction A2N1N2 <C A1N1N2, we have to leave out 
the corresponding A2 hopping term for consistency as well. 
This leads to the simplified steady-state equations 

= f(S l -N l )-W l N l -A l N l S 2 , 

= f(S 2 - A/2) - W 2 N 2 +A1N1S2 - 2AiNiN 2 , 

which yield an efficiency 

= 2MN1N2 = 2/A 1 S 1 S 2 (Vi +Ai5) 

nv fS S(V 1 +A l S2)[V2(V i +A 1 S 2 )+2fA l Siy 

(14) 

where V, = W, + /. We could now evaluate this at a tempera- 
ture right on the plateau. It will turn out, however, that we can 
make two more assumptions for this case. 
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FIG. 5: (Color online) Efficiency versus temperature for various fractions of deep sites. Left column AE = 250 K, middle AE = 750 K, right 
AE = 1500 K. Red circles: KMC longhop results (see Sec. IV). Red lines the numerical solution of rate equations with standard A; = a,/S 
(solid), blue lines contributions by reaction on the 1- and 2-sites (dashed), and by switching between the types (dot-dashed). Dotted green 
line the results of the plateau model, and green horizontal line the simple (16). Vertical green line at T eq . The first row shows the results for 
homogeneous systems of only standard or only deep sites, respectively. 



First, we also neglect desorption from the 2-sites, so V2 = /, S%/ S, we have V\ /a \ <C S2/S < 1 . This yields 
and using A\ =a\/S, Eq. (14) reduces to 2 

_ 2(s 1 /s)(s 2 /s)(i + v 1 / fll ) ^"sTfc+T (16) 

,p (Vi/ai+S 2 /S)(l + Vi/ai + Si/S)' '■ , , , , , , , 

v ' ' M ' ' ' which no longer depends on any energy scales, and which we 

Second, on the plateau and for a reasonable deep-site fraction find to be in excellent agreement with both KMC and full rate 
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equation results (Fig. 5): Whenever a plateau forms (i.e., if AE 
is large enough to separate the homogeneous-system peaks, 
and if there are enough deep wells if AE is fairly large), the 
above expression is valid. 

To check the validity of these approximations, we recall 
from Sec. IV B that the peak temperature r™ 2 * for standard- 
site parameters was found to always belong to the plateau. 
It is large enough not to lie on the low-temperature rise to 
the standard-site peak, yet minimal so as not to depend on 
the peak separation governed by AE. At T = r™ 2 ", V\ = 

Wi +f = 2f and V 2 = W 2 +f = /[(/ /v) AE/Ew i + 1]. Reason- 
ably, f/v 1, while the smallest interesting AE ~ Ew x — E ax , 
such that the ratio AE/E^ is not excessively smaller than 
unity. This justifies the approximation V 2 ~ /, immediately 
eliminating AE from the game, as suggested by Fig. 5. We 
now check the orderofVi/ai =2//ai =2(//v) (£M, i _£fl i )/£,v i. 
The exponent is about 0.22 for amorphous carbon, and with 
the corresponding standard flux we have V\/a\ ~ 6.3 x 10~ 5 
(cf. Sec. IV A). This is negligible compared to any interesting 
deep-well fraction 52/5, which completes the argument for 
Eq. (16). (We checked that this holds at least equally well for 
standard olivine parameters [27].) 

Knowing what terms can be neglected, this result is also eas- 
ily derived from further simplified rate equations. We rather 
provide an intuitive explanation. We consider the system in 
the steady state, so all particles entering the system also have 
to leave. They enter by impingement to any site, and leave 
only from 2-sites, by LH rejection or by reaction with an in- 
coming 1 -particle. Consequently, particles from 1 -sites arrive 
at a rate fS\ /S 2 at each 2-site. This implies a rate 2/5i /S 2 -N 2 
of particles to leave the system, as the reaction takes away two 
atoms. Alternatively, particles leave by LH rejection (rate- 
wise, this is merely a separate desorption process) at a rate 
f-N 2 . The efficiency is the fraction of impinging particles that 
react; in the steady state, this is just the rate at which particles 
leave due to reaction, normalized by the total rate to leave (by 
reaction or by LH rejection). This yields 

25 

^= 257^' 

which coincides with Eq. (16). We note that Eq. (17) can be 
rewritten as 

for 52/5 <C 1, which is precisely of the form of the empiri- 



cal relation (10). The coefficient inside the second bracket in 
Eq. (10) deviates from 1 /2 because it was obtained through a 
fit over the entire range of 52/5 € [0,1], whereas Eq. (18) is 
strictly valid only when 52/5 is small. 



VI. CONCLUSIONS 

We have studied diffusion-limited reactions of particles on 
a two-dimensional lattice which consists of shallow and deep 
sites, using KMC simulations and rate equations. In the case 
when the two types of sites are randomly mixed, we found 
that the temperature range in which the reaction is efficient 
dramatically broadens compared to a homogeneous system 
that includes only shallow or only deep sites. The rate equa- 
tions are found to provide a good description of the system 
and are in perfect agreement with the KMC results. We have 
also studied a system in which the deep sites are clustered to- 
gether. In this case the hopping between shallow and deep 
sites is suppressed. As a result, the recombination efficiency 
is dramatically reduced in comparison with the case in which 
the shallow and deep sites are randomly mixed. 

We expect that the qualitative features observed for the bi- 
nary distribution will hold for a broader class of models with 
different distributions of binding energies. The results pre- 
sented in this paper are also relevant in the context of molec- 
ular hydrogen formation in the interstellar medium. More 
specifically, high abundances of molecular hydrogen are ob- 
served in photon-dominated regions [37]. In these regions, 
the grain temperatures are too high to form molecular hydro- 
gen from weakly adsorbed hydrogen atoms. It was proposed 
that strong-binding sites in conjunction with the weak-binding 
sites enable the efficient formation of molecular hydrogen un- 
der these conditions [38, 39]. Our work provides a quantita- 
tive basis for this mechanism. 
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